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Abstract. We consider in this paper the multivariate regression problem, when the 
target regression matrix A is close to a low rank matrix. Our primary interest is in on 
the practical case where the variance of the noise is unknown. Our main contribution is 
to propose in this setting a criterion to select among a family of low rank estimators and 
prove a non-asymptotic oracle inequality for the resulting estimator. We also investigate 
the easier case where the variance of the noise is known and outline that the penalties 
appearing in our criterions are minimal (in some sense). These penalties involve the 
expected value of Ky-Fan norms of some random matrices. These quantities can be 
evaluated easily in practice and upper-bounds can be derived from recent results in 
random matrix theory. 



We build on ideas introduced in a recent paper of Bunea, She and Wegkamp [3 [8] for the 
multivariate regression problem 



where Y is a m x n matrix of response variables, X is a m x p matrix of predictors, A is 
p x n matrix of regression coefficients and E is a m x n random matrix with i.i.d. entries. 
We assume for simplicity that the entries Ei j are standard Gaussian, yet all the results 
can be extended to the case where the entries are sub-Gaussian. 

An important issue in multivariate regression is to estimate A or XA when the matrix 
A has a low rank or can be well approximated by a low rank matrix, see Izenman |14j . 
In this case, a small number of linear combinations of the predictors catches most of the 
non-random variation of the response Y. This framework arises in many applications, 
among which analysis of fMRI image data [11] , analysis of EEG data decoding [2\ , neural 
response modeling [6] or genomic data analysis [7]. 

When the variance a 2 is known, the strategy developed by Bunea et al. [7] for estimating A 
or XA is the following. Writing |.| for the Hilbert-Schmidt norm and A r for the minimizer 
of \\Y — XA\\ over the matrices A of rank at most r, the matrix XA is estimated by XA?, 



Date: September 2010, revision April 2011. 

2010 Mathematics Subject Classification. 62H99,60B20,62J05. 

Key words and phrases. Multivariate regression, random matrix, Ky-Fan norms, estimator selection. 



1. Introduction 




Y = XA + oE 
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where f minimizes the criterion 

(2) Crit ff2 (r) = \\Y-XA r \\ 2 + pen CT 2 (r)a 2 . 

Bunea et al. [7J considers a penalty pen CT 2 (r) linear in r and provides clean non-asymptotic 
bounds on H-X^ — XA|| 2 , on \\A? — A\\ 2 and on the probability that the estimated rank 
f coincides with the rank of A. 

Our main contribution is to propose and analyze a criterion to handle the case where a 2 
is unknown. Our theory requires no assumption on the design matrix X and applies in 
particular when the sample size m is smaller than the number of covariates p. We also 
exhibit a minimal sublinear penalty for the Criterion ^ for the case of known variance. 

Let us denote by q the rank of X and by G qxn a g x n random matrix with i.i.d. standard 
Gaussian entries. The penalties that we introduce involve the expected value of the Ky-Fan 
(2, r)-norm of the random matrix G qxn , namely 

r 

Sqxn(r) = E [|| Ggxn || (2,r)] : where ||Ggxn|| 2 2 ,r) = ^2 a k( G qxn) 

fe=l 

and where ak(G qxn ) stands for the k-th largest singular value of G qxn . The term S qxn {v*) 
can be evaluated by Monte Carlo and for q, n large enough an accurate approximation of 
Sqxn(r) is derived from the Marchenko-Pastur distribution, see Section 2. 

For the case of unknown variance, we prove a non-asymptotic oracle-like inequality for the 
criterion 

(3) Crit(r) = log(||Y - XA r \\ 2 ) + pen(r). 

with ^ 

pen(r) > - log ( 1 - K Sqxn ^ ) with K > I. 
\ nm — 1 / 

The latter constraint on the penalty is shown to be minimal (in some sense). In addition, 
we also consider the case where a 2 is known and show that the penalty pen(r) = S qxn {r) 2 
is minimal for the Criterion Q. 

The study of multivariate regression with rank constraints dates back to Anderson pQ 
and Izenman |13| . The question of rank selection has only been recently addressed by 
Anderson |T] in an asymptotic setting (with p fixed) and by Bunea et al. |H] in an non- 
asymptotic framework. We refer to the latter article for additional references. In parallel, 
a series of recent papers study the estimator A{ obtained by minimizing 

\\Y -XA\\ 2 + \Y^Ok{A) 

k 

see among others Yuan et al. [25] , Bach [3] , Neghaban and Wainwright [TU] , Lu et al. [T7] , 
Rohde and Tsybakov [20], Koltchinskii et al. [15J. Due to the penalty ^2 k crk(A), the 
estimator A{ has a small rank for A large enough and it is proven to have good statistical 
properties under some hypotheses on the design matrix X. We refer to Bunea et al. [8] for 
a detailed analysis of the similarities and the differences between A^ and their estimator. 
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Our paper is organized as follows. In the next section, we give a few results On Oqxny ) 
and on the estimator XA r . In Section 3, we analyze the case where the variance a 2 is 
known, which gives us a benchmark for the Section 4 where the case of unknown variance 
is tackled. In Section 5, we comment on the extension of the results to the case of sub- 
Gaussian errors and we outline that our theory provides a theoretically grounded criterion 
(in a non-asymptotic framework) to select the number r of components to be kept in a 
principal component analysis. Finally, we carry out an empirical study in Section 6 and 
prove the main results in Section 7. 

R-code. The estimation procedure described in sections 4 and 7 has been implemented in 
R. We provide the R-code (with a short notice) at the following URL : 
http : //www . cmap . polytechnique . f r/ ~giraud/ sof tware/KF . zip 

What is new here? The primary purpose of the first draft of the present paper |10j was 
to provide complements to the paper of Bunea et al. [?] in the two following directions: 

• to propose a selection criterion for the case of unknown variance, 

• to give some tighter results for Gaussian errors. 

During the reviewing process of the first draft of this paper, Bunea, She and Wegkamp 
wrote an augmented version of their paper [8] were they also investigate these two points. 
Let us comment briefly on the overlap between the results of these two simultaneous 
works |10l IB]. Let us start with the main contribution of our paper, which is to provide 
a selection criterion for the case of unknown variance. In Section 2.4 of [8], the authors 
propose and analyze a criterion to handle the case of unknown variance in the setting 
where the rank q of X is strictly smaller than the sample size m. In this favorable case, 
the variance a 2 can be conveniently estimated by 

o ||y-PY|| 2 

a = , with P the orthogonal projector onto the range of X, 

ran — qn 

which has the nice feature to be an unbiased estimator of a 2 independent of the collection 
of estimators {A r , r = 0, . . . , q}. Plugging this estimator a 2 in the Criterion S, Bunea et 
al. [8] proves a nice oracle bound. This approach no more applies in the general case where 
the rank of X can be as large as m, which is very likely to happen when the number p of 
covariates is larger than the sample size m. We provide in Section 4 an oracle inequality 
for the Criterion ^ with no restriction on the rank of X. 

Concerning the case of known variance : the final paper of Bunea et al. [8j proposes for 
Gaussian errors the penalty pen CT 2(r) = Kr{^fq + \fn) 2 with K > 1 which is close to 
ours for r <C min(q,n). For moderate to large r, we mention that our penalty ^ can be 
significantly smaller than r{y/q + \fn) 2 , see Figure [l] below. 

Notations. All along the paper, we write A* for the adjoint of the matrix A and cri(A) > 
02(^4) > . . . for its singular values ranked in a decreasing order. The Hilbert-Schmit norm 
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of A is denoted by \\A\\ = Tr(A* A) 1 / 2 and the Ky-Fan (2,r)-norm by 

/ r \V2 
W A W(?,r) = \J2 a k(A) 2 ) . 

Finally, for a random variable X, we write EpT] 2 for (E[A"]) 2 to avoid multiple parentheses. 

2. A FEW FACTS ON S qxn (r) AND X A r 

2.1. Bounds on S qxn (r). The expectation S qxn (r) = E [||G g xn||(2,r)] can De evaluated 
numerically by Monte Carlo with a few lines of R-code, see the Appendix. From a more 
theoretical point of view, we have the following bounds. 

Lemma 1. Assume that q < n. Then for any r < q, we have S qxn (r) 2 > r ( n — l/<z) and 
S q xn(r) 2 < min <j r (^fn + ^/q) 2 ,nq- ^ (y/n - Vk) 2 ,r + ^ i^/n + y 7 ? - k + 1 

L fe=r+l k=l ) 

When q > n the same result holds with q and n switched. In particular, for r = min(n, q), 
we have 

qn-l< S 2 xn (mm(n, q)) = E [||G gxn ||] 2 < qn. 

2 



The proof of the lemma is delayed to Section [7j The map r — > S qxn {r) and the up- 
per/lower bound of Lemma 1 are plotted in Figure 1 for q = 200 and n = 200 and 1000. 
We notice that the bound r — > S qxn {r) 2 < r(^/q + \/n) 2 looks sharp for small values of r, 
but it is quite loose for moderate to large values of r 

Finally, for large values of q and n, asymptotics formulaes for S qxn (r) can be useful. It 
is standard that when n,q go to infinity with q/n — > ft < 1, the empirical distribution 
of the eigenvalues of n~ l G qxn G* qxn converges almost surely to the the Marchenko-Pastur 
distribution [IB], which has a density on [(1 — \/]3) 2 , (1 + \//3) 2 ] given by 

As a consequence, when q and n go to infinity with q/n — > /3 < 1 and rjq — > a < 1, we 
have 

/•(l+v 7 ?) 2 

(4) S qxn (r) 2 ~nq xfp(x)dx, 

where x a is defined by 



/^(ar) dx = a. 

Since the role of q and n is symmetric, the same result holds when n/q — > (3 < 1 and 
r/n — > a < 1. This approximation Q can be evaluated efficiently (see the Appendix) 
and it turns to be a very accurate approximation of S qxn (r) for n, q large enough (say 
nq > 1000). 
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Figure 1. In bold red r — > S qxn (r) 2 , in solid black r — > r (y/n + yfq) 2 , 
in dashed blue the upper-bound of Lemma [TJ in dotted green the lower 
bound. Left: q = n = 200. Right: q = 200 and n = 1000. 

2.2. Computation of XA r . Next lemma provides a useful formula for XA r . 

Lemma 2. Write P for the projection matrix P = X(X*X) + X* , with (X*X) + the 
Moore-Penrose pseudo-inverse of X*X . Then, for any r < q we have XA r = (PY) r 
where (PY) r minimizes \\PY — B\\ 2 over the matrices B of rank at most r. 

As a consequence, writing PY = UT,V* for the singular value decomposition of PY, the 
matrix XA r is given by XA r = UY> r V* , where S r is obtained from S by setting (E r )j j = 
for i > r + 1 . 

Proof of Lemma g We note that \\PY-P(PY) r \\ 2 < \\PY-(PY) r \\ 2 and rank(P(PY) r ) < 
r, so P(PY) r = (PY) r . In particular, we have (PY) r = XA r , with A r = (X*X) + X*(PY) r . 
Since the rank of XA r is also at most r, we have 

||y-Xi r || 2 = \\Y - PY\\ 2 + \\PY - (PY) r \\ 2 

< \\Y - PY\\ 2 + \\PY - XA r \\ 2 = \\Y - XA r \\ 2 . 

Since the rank of A r is not larger than r, we then have A r = A r . 

3. The case of known variance 

In this section we revisit the results of Bunea et al. [TJ [8] for the case where a 2 is known. 
This analysis will give us a benchmark for the case of unknown variance. Next theorem 
states an oracle inequality for the selection Criterion ^ with penalty fulfilling pen^ (r) > 
KS a xn( r ) for K > 1. Later on, we will prove that the penalty pen CT 2(r) — S q xn( r ) is 
minimal in some sense. 
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Theorem 1. Assume that for some K > 1 we have 

(5) pen cr 2(r) > KS qxn (r) 2 for all r<mm(n,q). 
Then, when f is selected by minimizing the estimator A = Af satisfies 

(6) E \\XA-XA\\ 2 <c(K) mm {e^XA- XA r \\ 2 + pen a2 (r)a 2 + ct 2 | 
for some positive constant c(K) depending on K only. 



The risk bound (|6j) ensures that the risk of the estimator A is not larger (up to a constant) 
than the minimum over r of the sum of the risk of the estimator A r plus the penalty term 
pen CT 2(r)(7 2 . We will see below that this ensures that the estimator A is adaptive minimax. 

For r <C min(n, q), the penalty pen CT 2(r) = KS qxn {r) 2 is close to the penalty pen^. 2 = 
K(^/q + ^/n) 2 r proposed by Bunea et al. [8], but pen CT 2(r) can be significantly smaller than 
pen^. 2 (r) for moderate values of r, see Figure 1. Next proposition shows that choosing a 
penalty pen CT 2(r) = KS qxn (r) 2 with K < 1 can lead to a strong overfitting. 

Proposition 1. Assume that A = and that r is any minimizer of the Criterion |lp with 
pen CT 2(r) = KS qxn (r) 2 for some K < 1. Then, setting a = 1 — + K)/2 > we have 

A l-K nq-1 \ J/o e -« 2 ™«(n,g)/2 

r > x ; ^_ — > 1 - e ' 



4 (y/n + ^/g) 2 / 1 - e - " 2 max(n,q)/2 ' 

As a consequence, the risk bound |6p cannot hold when Condition |5p is replaced by 
pen CT 2(r) = KS qxn {r) 2 with K < 1. 

In the sense of Birge and Massart |5j, the Condition ^ is therefore minimal. 
Minimax adaptation. 

Fact 1. For any p S ]0, 1], there exists a constant c p > such that for any integers m, n,p 
larger than 2, any positive integer q less than min(m,p) and any design matrix X fulfilling 

(7) a q (X) > pai(X), where q = rank(X), 

we have 



inf sup E 

A A:rank(A)<r 



\XA-XA\\ 2 



> c p (q + n)ro~ 2 , for all r < min(n, q). 



When p < m and q = p, this minimax bound follows directly from Theorem 5 in Rohde 
and Tsybakov [20] as noticed by Bunea et al, see [8] Section 2.3 Remark (ii) for a slightly 
different statment of this bound. We refer to Section 7/7 for a proof of the general case 
(with possibly q < p and/or p > m). 

If we choose pen CT 2(r) = KS qxn (r) 2 for some K > 1, we have pen cr 2(r) < 2Kr(q + n) 
according to Lemma 1. The risk bound (JgJ) then ensures that our estimator A is adaptive 
minimax (as is the estimator proposed by Bunea et al.). 
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4. The case of unknown variance 



We present now our main result which provides a selection criterion for the case where 
the variance a 2 is unknown. For a given r max < min(rt, q), we propose to select f 6 
{1, . . . , r max } by minimizing over {1, . . . , r max } Criterion namely 

Crit(r) = log(||Y - XA r \\ 2 ) + pen(r). 

We note that the Criterion ^ is equivalent to the criterion 

2 f-, , pen'(rY 



(8) 



Crit'(r) = ||y-X^ r || 2 1 



nm 



with pen'(r) = nm(e pen ^ — 1). This last criterion bears some similitude with the Crite- 
rion ([2]). Indeed, the Criterion ^ can be written as 



\Y - XA r \\ 2 + pen'(r)a 2 , 



with a 2 = \\Y — X A r \\ 2 / (nm) , which is the maximum likelihood estimator of a 2 associated 
to A r . To facilitate comparisons with the case of known variance, we will work henceforth 
with the Criterion Q . Next theorem provides an upper bound for the risk of the estimator 
XA f . 



Theorem 2. Assume that for some K > 1 we have both 

■ ,.nax) 2 + 1 < nm 

\2 



(9) 
(10) 



KSqxn( r n 



and pen'(r) > 



KSq X n( r 



1 



1 

nm 



(1 + KS, 



qxn 



r) 



5 



Then, when r is selected by minimizing 
satisfies 



for r < r max . 
over {1, . . . , r max }, the estimator A = A? 



(11) E 



\XA-XA\ 



\XA r -XA\\ 2 
for some constant c(K) > depending only on K . 



< c{K) min { E 



I + ) + (pen' (r) + l)a 2 



nm 



Let us compare Theorem [2] with Theorem [T] The two main differences lie in Condition (10) 
and in the form of the risk bound ( 11 ) . Condition ( 10 ) is more stringent than Condition |5^. 
More precisely, when r is small compared to q and n, both conditions are close, but when 
r is of a size comparable to q or n, Condition (10) is much stronger than ([5]). In the case 
where m = q, it even enforces a blow up of the penalty pen'(r) when r tends to min(n, m). 
This blow up is actually necessary to avoid overfitting since, in this case, the residual 
sum of squares \\Y — XA r \\ 2 tends to when r increases. The second major difference 
between Theorem [2] and Theorem [I| lie s in the multiplicative factor ( 1 + pen' (r)/nm) in the 
right-hand side of the risk bound (11). Due to this term, the bound (11) is not (strictly 



s 
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speaking) an oracle bound. To obtain an oracle bound, we have to add a condition on r max 
to ensure that pen'(r) < Cnm for all r < r max . Next corollary provides such a condition. 

Corollary 1. Assume that 

nm — 1 

(12) r miD[ < a — t — — for some < a < 1, 



and set 



pen(r) = — log(l — KS q x n {f) 2 1 '{nm — 1)) for some K > 1. 

Then, there exists c(K, a) > such that, when f is selected by minimizing 
{1, . . . ,r max }, we have the oracle inequality 



over 



E 



\XA-XA\\ 2 



< c(K,a) min [l 



\XA r - XA\\ 



In particular, the estimator A is adaptive mininicLX up to th.6 rank r m ax 

specified by al2«. 



In the worst case where m = q, Condition (12) requires that r max remains smaller than a 



fraction of min(n, q). In the more favorable case where m is larger than Aq, Condition ( 12 ) 
can be met with r max = min(g, n) for suitable choices of K and a. 



n=q=200 



n=1000, q=200 





Figure 2. In dotted green pen(r) = — log(l — S qxn (r) 2 /(nq — 1)), in 
solid black pen(r) = r (yjn -\- ^/q) 2 /{nq) : in dashed red pen' (r)/(nq) = 
S q xn(r) 2 /{nq - 1 - S qxn {r) 2 ). Left : q = n = 200. Right : q = 200 and 
n = 1000. 
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Let us discuss now in more details the Conditions (| 9j) and ( 10 ) of Theorem [2] We have 



S q xn{f) < r(y/n + y/q) so the Conditions (191) and (10) are satisfied as soon 



as 



nm-1 , Kr{^h+^/q) 
—. —k and pen (r) > -, — — — — , for r < r„ 



In terms of the Criterion ([3]), the Condition (10) reads 

pen(r) > - log(l - KS qxn (r) 2 / {nm - 1)). 

When pen(r) is defined by taking equality in the above inequality, we have pen(r) ~ 
Kr(^Jn + yfq) 2 / {nm) for small values of r, see Figure [2J 

Finally, next proposition, shows that the Condition ( |10[ ) on pen'(r) is necessary to avoid 
overfitting. 

Proposition 2. Assume that A = and that f is any minimizer of Criterion |5p over 
{1, . . . , min(n, q) — 1} with 



/ 10 \ // \ KS q xn{ r ) 2 , - 

(13) pen (r) = ^— — for some K < 1. 

Then, setting a = (1 — K)/A > we have 

l-K nq-l \_ _ ^2 /0 e -« 2 max(n >f/ )/2 



(Vn + y/q) 2 J~ 1 _ e -a 2 max(n,g)/2 ' 

As in Proposition [TJ a consequence of Proposition [2] is that Theorem [2] cannot hold with 



Condition (10) replaced by (13). Condition (10) is then minimal in the sense of Birge and 
Massart [5]. 

5. Comments and extensions 

5.1. Link with PC A. In the case where X is the identity matrix, namely Y = A + E, 
Principal Component Analysis (PCA) is a popular technique to estimate A. The matrix 
A is estimated by projecting the data Y on the r first principal components, the number 
r of components being chosen according to empirical or asymptotical criterions. 

It turns out that the projection of the data Y on the r first principal components coincides 
with the estimator A r . The criterion (J3J) then provides a theoretically grounded way to 
select the number r of components. Theorem [2] ensures that the risk of the final estimate 



Af nearly achieves the minimum over r of the risks E 



\A r 



5.2. Sub-Gaussian errors. We have considered for simplicity the case of Gaussian er- 
rors, but the results can be extended to the case where the entries E% j are i.i.d sub- 
Gaussian. In this case, the matrix PE will play the role of the matrix G q xn in the Gauss- 
ian case. More precisely, combining recent results of Rudelson and Vershynin |21j and 
Bunea et al. [7'| on sub-Gaussian random matrices, with concentration inequality for sub- 
Gaussian random variables M enables to prove an analog of Lemma [l] for E [\\PE k ^ 



10 CHRISTOPHE GIRAUD 

(with different constants). Then, the proof of Theorem [I] and Theorem [2] can be easily 
adapted, replacing the Condition ([5]) by 

pen(r) > KK [\\PE\\(2 >r )\ 2 , for r < mm(q,n), 



and the Conditions and |lu| by ATE [11-^11 (2,rw)] 2 < E[||£||] 2 and 

1 2 



KE l\\PE\\( 2 ,r)\ f < 

pen (r) > k -, tor r < r n 

l-KE[\\PE\\ {2 , r) ] 2 /E[\\E\\} 2 

Analogs of Proposition 1 and 2 also hold with different constants. 



5.3. Selecting among arbitrary estimators. Our theory provides a procedure to select 
among the family of estimators {A r , r < r max }. It turns out that it can be extended to 
arbitrary (finite) families of estimators {A\, A G A} such as the nuclear norm penalized 
estimator family {A^ , A G A}. The most straightforward way is to replace everywhere 
A r by A\ and pen(r) by pen(A), with pen(A) = pen(rank(^A))- In the spirit of Baraud et 
al. [1], we may also consider more refined criterions such as 

Crit Q (A) = min 1{\\Y - XA x , r \\ 2 + \\XA X - XA X J 2 ) ( 1 + P<, "' ( '' ' 



r<r max [ y nm 

where a > and A\^ r minimizes \\B — A\\\ over the matrices B of rank at most r. Analogs 
of Theorem [2] can be derived for such criterions, but we will not pursue in that direction. 

6. Numerical experiments 

We perform numerical experiments on synthetic data in two different settings. In the first 
experiment, we consider a favorable setting where the sample size m is large compared to 
the number p of covariables. In the second experiment, we consider a more challenging 
setting where the sample size m is small compared to p. The objectives of our experiments 
are mainly: 

• to give an example of implementation of our procedure, 

• to demonstrate that it can handle high-dimensional settings. 



Simulation setting. Our experiments are inspired by those of Bunea et al. [3 [8] , the 
main difference is that we work in higher dimensions. The simulation design is the follow- 
ing. The rows of the matrix X are drawn independently according to a centered Gaussian 
distribution with covariance matrix Ejj = p'* - ^', p > 0. For a positive b, the matrix A is 
given by A = bB pxr B rxn , where the entries of the B matrices are i.i.d. standard Gaussian. 
For r < mm(n,p), the rank of the matrix A is then r with probability one and the rank 
of X is min(m,p) a.s. 
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Experiment 1: in the first experiment, we consider a case where the sample size m = 
400 is large compared to the number p = 100 of covariables and n = 100. The other 
parameters are r = 40, p varies in {0.1,0.5,0.9} and b varies in {0.025,0.05,0.075,0.1}. 
This experiment is actually the same as the Experiment 1 in [8], except that we have 
multiplied m, p, n, r by four and adjusted the values of b. 

Experiment 2: the second experiment is much more challenging since the sample size 
m = 100 is small compared to the number p = 500 of covariables and n = 500. Fur- 
thermore, the rank q of X equals m, which is the least-favorable case for estimating the 
variance. For the other parameters, we set r = 20, p varies in {0.1,0.5,0.9} and b varies 
in {0.005,0.01,0.015,0.02}. 

Estimators. For K > 1, we write KF[_RT] for the estimator A? with f selected by the 
Criterion ([8]) with 

i, \ KS qxn (r) 2 

Pen(r) " 1-^(1 + KS qxn {rY) 
(the notation KF refers to the Ky-Fan norms involved in S qxn {r)). 

For A > 0, we write RSC[A] for the estimator A? with f selected by the criterion introduced 
by Bunea, She and Wegkamp [7] 

Crit A (r) = \\Y - XA r f + A(n + rank(X))r. 

Bunea et al. [8J proposes to use A = Ka 2 with K > 2 and 

\\Y — PY\\ 2 

a 2 = — -, with P the projector onto the range of X. 

mn — nrank(X) 

We denote by RSCI[iT] the resulting estimator RSC[K<r 2 ]. 

Both procedures KF and RSCI depend on a tuning parameter K. There is no reason for the 
existence of a universal "optimal" constant K. Nevertheless, Birge and Massart [5] suggest 
to penalize by twice the minimal penalty, which corresponds to the choice K = 2 for KF. 
The value K = 2 is also the value recommended by Bunea et al. [8\ Section 4 for the RSCI 
(see the "adaptive tuning parameter" fJ, a dapt)- Another classical approach for choosing a 
tuning parameter is Cross- Validation : for example, K can be selected among a small 
grid of values between 1 and 3 by V-Fold CV. We emphasize that there is no theoretical 
justification that Cross- Validation will choose the best value for K. Nevertheless, for 
each value K in the grid, the estimators KF[if] and RSCIf-fT] fulfills an oracle inequality 
with large probability, so the estimators with K chosen by CV will also fulfills an oracle 
inequality with large probability (as long as the size of the grid remains small). We will 
write YF[K = CV] and RSCI[K = CV] for the estimators KF and RSCI with K selected 
by 10-fold Cross- Validation. 

Finally, in Experiment 2 the estimator RSCI cannot be implemented since rank(X) = m. 
Yet, it is still possible to implement the procedure RSC[A] and select A > among a 
large grid of values by 10-fold Cross- Validation, even if in this case there is no theoretical 
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control on the risk of the resulting estimator RSC[A = CV]. We will use this estimator as 
a benchmark in Experiment 2. 

Results. The results of the first experiment are reported in Figure [3] and those of the 
second experiment in Figure |4j The boxplots of the first line compare the performances of 
estimators KF and RSCI to that of the estimator XA r that we would use if we knew that 
the rank of A is r. The boxplots give for each value of p the distribution of the ratios 

, , \\xa-xa\\ 2 

(14) J! ^J— A 10 

\\XA-XA T \\ 2 

for A given by KF[K = 2], RSCI[# = 2], KF[K = CV] and RSCI[if = CV] in the first 
experiment and by KF[if = 2], KF[K = CV] , and RSC[A = Cy] in the second experiment. 



The ratios (14) are truncated to 10 for a better visualization. Finally, we plot in the 



second line the mean estimated rank E[r] for each estimator and each value of b and p. 
Experiment 1 (large sample size). All estimators VF[K = 2], RSCl[K = 2], ¥F[K = CV] 



and RSCI[-fr = CV] perform very similarly and almost all the ratios (14) are equal to 1. 



Experiment 2 (small sample size). The estimator KF[.K" = CV] has global good perfor- 



mances, with a median ratio (14) around 1, but the ratio (14) can be as high as 5 in 
some examples for p = 0.9. In contrast, the estimator KF[iT = 2] is very stable but it 
has a median value significantly above the other methods. Finally, the performances of 
the estimator RSC[A = CV] are very contrasted. For small correlation (p = 0.1), its per- 
formances are similar to that of KF[_?T = CV]. For p = 0.5 or p = 0.9, it has very good 
performances most of the time (similar to KF[A" = CV]) but it completely fails on a small 



fraction of examples. For example, for p = 0.9, it has a ratio (14) smaller than 7 in 80% of 
the examples (with a median value close to 1), but in 20% of the examples, it completely 
fails and the ratio \\XA - XA\\ 2 /\\XA - XA r \\ 2 for RSC[A = CV] can be has high as 10 13 
(these values do not appear in Figure [4] since we have truncated the ratio ( 14 ) to 10 to 
avoid a complete shrinkage of the boxplots). We recall, that there exists no risk bound 
for the estimator RSC[A = CV], so these results are not in contradiction with any theory. 

Finally, we emphasize that no conclusion should be drawn from these two single experi- 
ments about the superiority of one procedure over the others. 



7. Proof of the mains results 

7.1. Proof of Lemma [T| For notational simplicity we write G = G qX n- The case r = 1 
follows from Slepian's Lemma, see Davidson and Szarek [9] Chapter 8. For r > 1, we note 
that 

E [l|G|| (2 , r) ] 2 < min jrE [||G|| (2>1) ] 2 , £ E[a 2 k (G)] j . 
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The first upper bound S qxn (r) 2 < r(yjn + ^/q) 2 follows. For the second upper bound, we 
note that 

]TE[^(G)]<E[||G|| 2 ]- jz E[a k (G)) 2 . 

fc=l k=r+l 

The interlacing inequalities [12] ensure that a k (G' k ) < cr k (G) where G' k is the matrix made 
of the k first rows of G. The second bound then follows from E[cr/%(G" fc )] > y/n — \/k, 
see [9]. 

Let us turn to the third bound. The map G —> o~ k (G) is 1-Lipschitz so, writing M k for 
the median of a k , the concentration inequality for Gaussian random variables ensures that 
(Mfc — a k (G)) + < £ + and (a k (G) — M&)+ < where £+ and are the positive part of 
two standard Gaussian random variables. As a consequence we have 

E[a 2 (G)] - E[a k (G)] 2 + (M k - E[a k (G)}) 2 = E [(a k (G) - M k )%] + E [{M k - a k {G))\] 

< E[c >2 ]+E[e + ] = i, 

and thus E[a 2 (G)} < E[a k {G)} 2 + 1. 

Furthermore, the interlacing inequalities [T2] ensure that a k (G) < a\{G'_ We can 
then bound E[cjjt(G)] by 

E[a k (G)} < sfn + y/q - k + 1 
which leads to the last upper bound. 

For the lower bound, we start from HGH^r) — ll ( -^ll 2r /'? ( sum of a decreasing sequence) 
and use again the Gaussian concentration inequality to get 

E[||G|| 2 ]-l = n<7-l<E[||G||] 2 
and concludes that r(nq — l)/q < E [||G||( 2 ,r)] 2 = <S qxn (r) 2 . 



7.2. A technical lemma. Next lemma provides a control of the size of the scalar product 
< E,XA k — XA r > which will be useful for the proofs of Theorem [l] and Theorem [2] 

Lemma 3. Fix r < min(ra, q) and write A r for the best approximation of A with rank at 
most r. Then, there exists a random variable U r such that E(U r ) < rmin(n, q) and for 
any rj > and all k < min(n, q) 

~ 1 ~ l + l/r/ 

(15) 2a I < E, XA k - XA r > I < \\XA k - XA\\ 2 + . '.' \\XA - XA r \\ 2 

1 + 7? (l + r?) 2 

+ (1 + t/) 2 (1 + l/v)a 2 U r + (1 + nfa 2 \\PE\\ 2 m 



where P = X(X* X) + X* is as in Lemma^ 
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Iterating twice the inequality lab < a 2 /c + cb 2 gives 
2a\ < E, XA k - XA r > | 

< -^-\\XA k - XAf + ^> 2 \\XA - XA r f + (1 + vfa 2<E,XA k -XA r > 2 
~1 + V U " (1 + ??) 2 " " V /; \\XA k -XA r \\ 2 

We write XA r = UT r V* for the singular value decomposition of XA r , with the convention 
that the diagonal entries of T r are decreasing. Since the rank of XA r is upper bounded by 
the rank of A r , the m x n diagonal matrix T r has at most r non zeros elements. Assume 
first that n < q. Denoting by I r the m x m diagonal matrix with (I r )a = 1 if i < r and 
(Ir)i,i = if i > r and writing I_ r = I — I r and B k = U* X A k V , we have 

< S, - Ii r > 2 < U*PEV, B k - T r > 2 



\XA k -XA r \\ 2 ||s fc -r r || 2 

< V"' PEY.IABi, - I.Y) > + < U*PEV,I- r B k E' 

2 . II r R, 112 



|I r (S fc -r r )|| 2 + ||/_ rJ B fc | 



^ _! < C/*PW,/ r (5 fc -r r ) > 2 , , .<U*PEV,I- r B k >< 

< {l + V ) ~ 1" (1 + T])- 



\E(B k -T r )\\ 2 \\I~rB k \\ 2 
The first term is upper bounded by 

< U*PEV,I r (B k -T r ) > 2 „ „ ll9 
1 - < \I r U*PEV 2 = U r 

\\I r (B k -T r )\\ 2 

and the expected value of the right-hand side fulfills 

E(U r ) = n\\I r U*PU\\ 2 = n\\U*PUI r \\ 2 < nr. 

Since the rank of I- r B k is at most k, the second term can be bounded by 

<U*PEV,I^ r B k > 2 <U*PEV,B> 2 _ 2 _ 2 
1, r s n2 sup iirP _ W u FhV W{2,k) ~ \\ Ftj W%ky 

\\I- r B k \\ 2 rank(B)<fc \\&\\ 



Putting pieces together gives (15) for n < q. The case n > q can be treated in the same 
way, starting from 

B k — T r = (B k — T r )I r + B k I- r 
with I r and I- r two n x n diagonal matrices defined as above. 

7.3. Proof of Theorem [Tj The inequality Crit (T 2(f) < Crit cr 2(r) gives 

(16) \\XA-XA\\ 2 < \\XA r - XA\\ 2 + pen a2 {r)a 2 + 2a < E,XA-XA r > -pen^f^ 2 . 
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Combining this inequality with Inequality (15) of Lemma|3|with 7] = ((l+ET)/2) 1 / 3 — 1 > 0, 
we obtain 



—9— \\XA-XA\\ 2 < \\XA r - XA\\ 2 + 2 1 + ^ \\XA - XA r f + 2pen ff2 (r)a 2 



+2(1 + 7?) 2 (1 + »T V^r + ^tl^p^^ - pen ff2 (r)a 2 

min(n,<?) 



fe=i ^ ' - 



The map £7 — )• ||P£?[|/ 2j fe) is 1-Lipschitz and convex, so there exists a standard Gaussian 
random variable £ such that [|P£?[L 2) fc) < ^[||-f-^l|(2,Jfc)] + £+ an< ^ then 

^±1 ||PP|| 2 2 , fc) - pen(fc)) < + 2t + E[\\PE\\ {2jk) } - ^ziE[||PP|| (2>fc) ] 2 ) 



E 



< Cl (K)eM-C2(K)E[\\PE\\ m } 



2\ 



Since ||PP||(2,fc) is distributed as ||G g xn||(2,fc)> Lemmajljgives that E[||PE?||(- 2)fe )] 2 > fcmax(n, q)- 
1 and the series 

mm(n,q) / 

\PEf m - V eu{k) 



fe=i v 



can be upper-bounded by Cl (K)e C2 ^ (l - e~ C2 ^) ^^WHm), Finally, E[U r ] < 



r min(n, g) is bounded by l + pen(r) and ||Xt4 — XAj-H 2 is smaller than E 
so there exists some constant c(K) > such that ^ holds. 



\XA-XAr\\ 2 



7.4. Proof of Theorem [2j To simplify the formulaes, we will note pen(r) = pen' (r)/(nm). 
The inequality Crit'(f) < Crit (r) gives 

\\XA- XA\\ 2 (l + pen(f)) 

< \\Y - XA r \\ 2 - a 2 (l + pen(r))||P|| 2 + pen(r)||y - XA r \\ 2 + pen(r)||P|| V 
+2(1 + pen(f))o- < E,XA-XA > -pen(f) ||P|| 2 o- 2 

< (2a < E,XA r - XA r > -pen(r)a 2 ||P|| 2 ) + + (1 + 2peH(r))||XA - XA r \\ 2 + 3pen(r) ||P|| 2 

+(1 + pen(r)) ( 2a < E, XA - XA r > P en ( r ) ; \\E\\ 2 a 2 ) + 2crpen(f) < E, XA r - XA > . 

V l + pen(r) J + 

Dividing both side by 1 + pen(f), we obtain 

\\XA-XA\\ 2 < (l+2pen(r))||X^-Xl r || 2 +3pen(r)||P|| 2 +2a| < E,XA-XA r > \+A r +A f . 
where 

A, ( 2a\ < E,XA k - XA r > \ - - ggg- ||P|| V 

1 + pen(fc) 
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We first note that E[||J5|R = nm and 2aE[\ < E,XA - XA r > |] < a 2 + \\XA - XA r \\ 2 . 
Then, combining Lemma pi with rj = (K 1 ^ — 1) and the following lemma with 5 = n gives 



E 



\XA-XA\\ 2 



< c(K) IE 



(1 + pen(r)) + (1 + nmpen(r))cr 2 ) , 



\XA-XA r \\ 2 
for some c{K) > 0. 

Lemma 4. Write P for the projection matrix P = X(X*X) + X* , with {X*X) + the 
Moore-Penrose pseudo-inverse of X*X. For any 5 > and r < min(n,q) such that 
(1 + 6)E [\\PE\\(2 r )] < Vnm - 1, we have 
(17) 



E 



PE\\fa r) - (1 + SfE [\\PE\\ {2jr) \ 2 \\E\\ 2 /(nm - 1))+ < 4 (1 + 1/5) e<^ 4 e 



= <5 2 /4 — 8 2 r max(n,<j)/4 



As a consequence, we have 



E 



sup (\\PE\\ 2 ^ r) - (1 + 5fS qxn (r) 2 \\E\\ 2 /(nm - 1] 

T^max 



2 1 A e 



-<5 2 max(n,g)/4 



< 4 (1 + 1/5) e d /4 — T2 7 -^- n . 

Proof of the Lemma. 

Writing t = (1 + «J)E[||Pf7||( 2ir )]/E[||£;||] < 1, the map E -> ||P-E||( 2 ,r) - is \/2- 

Lipschitz. Gaussian concentration inequality then ensures that 

\\PE\\ %r) < t\\E\\+E[\\PE\\ i2tr) -t\\E\\] + 2y/Z 
< i||^|| + (2^-^E[||PE|| (2ir) ]) + , 
with £ a standard exponential random variable. We then get that 

\\PE\\\ 2tr) < (1 + 5)t 2 \\E\\ 2 + 4(1 + 1/5) (v/l - 5E[\\PE\\ {2 , r) ]/2^ * 



and 



E 



PE\\ 2 {%r) - (l + 5)t z \\E\\' z ) + < 4(1 + 1/«5)E 



j.2i| cii|2 



^-5E[\\PE\\ {2; , 



< A(l + l/5)e- s2 ^ PE k^?/\ 



The bound ( 17) then follows from E[||Pi?||(2 jr )] 2 > r max(n, q) — 1 and E[||.E||] 2 > nm — 1. 
□ 

7.5. Proof Proposition [ll For simplicity we consider first the case where m = q. With 
no loss of generality, we can also assume that o~ 2 = 1. We set 

min(n,m) 

= {\\E\\ > (1 - a)E[\\E\\]} f| {||£|| (2>r) < (1 + a )E[\\E\\ (2ir) )} . 
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According to the Gaussian concentration inequality we have 

min(n,m) 

P(«o) > 1- E e- a2 ™W)l 2 / 2 
r=l 

min(n,m) 
a 2 /2 \ g-« 2 rmax(n,m)/2 



> i_ e «Va 



r=l 



where the last bound follows from Lemma [lj Furthermore, Lemma [i] gives that xA r = 
Y r {= E r ), where Y r is the matrix M minimizing \\Y — M\\ 2 over the matrices of rank at 
most r. As a consequence, writing m* = min(n,m), we have on £Iq 

Crit CT2 (m*) - Oritur) = KE^EW^] 2 - (\\E\\ 2 (2, m *) ~ \\ E \\ 2 (2,r)) ~ KE[\\E\\(2,r)] 2 

< ((i + a f _ K ) E[||i?|| (2ir) ] 2 - ((1 - a) 2 - K) EIII^H^)] 2 

< 2E[||^|| (2ir) ] 2 -i^E[|| J E|| (2im * ) ] 2 

< 2r(Jn + ^m) (nm — 1). 



We then conclude that on Qn we have r > x j_ m 

u — 4 ( v / n+- v /m)^ 

Let r* be the smaller integer larger than 1=^ x " m -j_ . Since IIXA— XA\\ 2 = \\E\\% ~\, 
we have 



E 



|Xl-Xy4|| 2 



> E ||i5/|| 2 2 jr *-)lf> r * 

> (1 - a) 2 S mxn (r*) 2 P ({r > r*} n {||£||( 2 ,r«) > (1 - a)5 m xn(r*)}) . 



Combining the analysis above with Gaussian concentration inequality for ||£J|b 2 ,r*)) we 
have 



g -a 2 max(n,m)/2 



P ({r > r*} n > (1 - a) W*)}) > 1 - 2e« / 2 f _ e _ Q2 max(n , m)/2 • 

We finally obtain the lower bound on the risk 



E 



\XA- XA\ 



> (1 - a)V(max(n,m) - 1) 1 - 2e a /2 



2/0 e 



-a 2 max(n,m)/2 



I g-a 2 max(n,m)/2 ; ' 



which is not compatible with the upper bound c(K) that we would have if ^ were also 
true with K < 1. 

When q < m, we start from ||Y - xl r || 2 = \\Y - PY\\ 2 + \\PY - XA r \\ 2 with P = 
X{X* X) + X* and follow the same lines, replacing everywhere E by PE and m by q. 
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7.6. Proof of Proposition [2j As in the proof of Proposition [TJ we restrict for simplicity 
to the case where a 2 = 1 and q = m, the general case being treated similarly. We write 
pen(r) = pen' (r)/(nm) and for any integer r* E [min(ra, m)/2, min(n, m) — 1], we set 

min(n,m) 

a* = {\\E\\(2,r*) > (1 - a)E[||£|| (2)r , } ]} f| {||£7[| (2 , r) < (1 + a)E[\\E\\ M )} . 

r=l 

According to the Gaussian concentration inequality we have 

min(n,m) 

P(0*) > 1-2 £ e ^[||s|| (2 ,,)] 2 /2 

r=l 

min(n,m) 

^ 2e a2 ^ ^ ^ g— a 2 r max(n,m)/2 

r=l 

where the last bound follows from Lemma [T] For any r < r*, we have on f2* 

Crit'(r*) - Crit'(r) = ||£7|| 2 (peH(r*) - pen(r)) + ||£|| 2 2r) (l + pen(r)) - \\E\\ 2 {2r) (l + pen(r*)) 

< (1 + a) 2 (E[||£||] 2 (plm(r*) - pen(r)) + E[||^|| (v) ] 2 (l + pen(r))(l + a) 2 
-Elll^^a + pen^^a-a) 2 . 

Since E[||£||] 2 < nm = KE[||£|| (2)r) ] 2 (l + pen(r))/pen(r), we have 

CritV) - Crit'(r) < (1 + a) 2 (l - K)(l + pn(r))E[||£|| (2jr) ] 2 

_((! _ a) 2 - (1 + a) 2 K){l + pm(r*))E[||£|| (2)r .)] 2 

< (1 + a) 2 (l - + pen(r*)) [E[||£|| (2>r) ] 2 - (1 - (1 + a)- 2 )E[|| J B|| (2>r , ) ] 2 ] 

To conclude, we note that E[||£||( 2ir )] 2 < r( v / n + \fm) 2 , E[||.E||( 2)r *)] 2 > (nm - l)/2 and 
1 — (1 + a)~ 2 > a, so the term in the bracket is smaller than 

r[y/n + ym) (nm—1) 

8 

which is negative when r < x ""' 1 



7.7. Minimax rate : proof of Fact 1. Let X = UY.V* be a SVD decomposition of 
X, with the diagonal elements of E ranked in decreasing order. Write U q and V q for the 
matrices derived from U and V by keeping the g-first columns, and Eg for q x q upper-left 
block of E (with notations as in R, U q = U[ , 1 : q], V q = V[ , 1 : q] and Eg = E[l : q, 1 : <?]). 
We have X = U q Y, q V* and 

y = ZB + aE, with Z = U q Y, q G R mxg and B = V* A € R" xn . 
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Let A be an arbitrary estimator of A and set B = V* A. Write Z* for the ith row of Z 
and {ei, . . . , e n } for the canonical basis of R n . According to ([7]), the map 

B^C(B) = [<Z ie * a ,B> / y 



a = 1, . . . , n 

fulfills the Restricted Isometry condition RI(r, v) of Rohde and Tsybakov [20] for all r < 
min(n, q) with 

2 2mn a x 1 ~ P 2 / 1 
^ = — 7^7^ and o = ^ < 1. 

CTl(X) 2 +f7 ? (X) 2 l + p 2 

Theorem 2.5 in [20j (with a = 1/10 and A = +oo) then ensures that there exists some 
constant c p > depending only on p such that 

inf sup E[||Zi? — ZB\\ 2 ] > 2c p (q + n)ra 2 , for all r < min(n, g). 

B B:vank(B)<r 

Let 5' be such that E[\\ZB - ZB'\\ 2 ] > c p (q + n)ro 2 and rank( J B') < r. The matrix 
A' = V q B' £ RP xn fulfills rank(^') < r and 



E 



\XA-XA'\\ 2 



E 



> c p (q + n)ra 2 . 



In conclusion, for any X fulfilling (I7j), any estimator A and any r < min(n, q), we have 



sup E 

A:rank(A)<r 



\XA-XA\\ 



> c p (q + n)ra 2 . 



References 



[i] 

[2] 

[3] 

[4] 
[5] 

[6] 

[7] 

[8] 

[9] 

[10] 
[11] 



T.W. Anderson. Estimating linear restrictions on regression coefficients for multivariate normal dis- 
tribution. Annals of Mathematical Statistics 22 (1951), 327-351. 

C. W. Anderson, E. A. Stolz, and S. Shamsunder. Multivariate autoregressive models for classication 
of spontaneous electroencephalogram during mental tasks. IEEE Trans, on bio-medical engineering, 
45 no 3 (1998), 277-286. 

F. Bach. Consistency of trace norm minimization, Journal of Machine Learning Research, 9 (2008), 

1019-1048. 

Y. Baraud, C. Giraud and S. Huet. Estimator selection in the Gaussian setting. |arXiv:1007. 2096 iT 
L. Birge and P. Massart. Minimal penalties for Gaussian model selection. Probability Theory and 
Related Fields, 138 no 1-2 (2007), 33-73. 

E. N. Brown, R. E. Kass, and P. P. Mitra. Multiple neural spike train data analysis: state-of-the-art 
and future challenges. Nature Neuroscience, 7 no 5 (2004), 456-461. 

F. Bunea, Y. She and M. Wegkamp. Adaptive rank Penalized Estimators in Multivariate Regression. 
|arXiv: 1004.2995^ 1 (2010) 

F. Bunea, Y. She and M. Wegkamp. Optimal selection of reduced rank estimation of high-dimensional 
matrices. To appear in the Annals of Statist. 

Davidson and Szarek. Handbook of the Geometry of Banach Spaces. North-Holland Publishing Co., 
Amsterdam, 2001. 

C. Giraud. Low rank multivariate regression. |arXiv:1009.5165V l (Sept. 2010) 



L. Harrison, W. D. Penny, and K. Friston. Multivariate autoregressive modeling of fmri time series. 
Neurolmage, 19 (2004), 1477-1491 



20 



CHRISTOPHE GIRAUD 



[12] R.A. Horn and C.R. Johnson. Topics in Matrix Analysis. Cambridge University Press, Cambridge, 
1994. 

[13] A.J. Izenman. Reduced-rank regression for the multivariate linear model. Journal of Multivariate 

analysis 5 (1975), 248-262. 
[14] A.J. Izenman. Modern Multivariate Statistical Techniques: Regression, Classification, and Manifold 

Learning. Springer, New York, 2008. 
[15] V. Koltchinskii, A. Tsybakov and K. Lounici. Nuclear norm penalization and optimal rates for noisy 

low rank matrix completion. arXiv:1011.6256v2 
[16] M. Ledoux. The concentration of measure phenomenon. Mathematical Surveys and Monographs, 89. 

American Mathematical Society, Providence, 2001. 
[17] Z. Lu, R. Monteiro and M. Yuan. Convex optimization methods for dimension reduction and coefficient 

estimation in multivariate linear regression. Mathematical Programming (to appear). 
[18] V. A. Marchenko, L. A. Pastur. Distribution of eigenvalues for some sets of random matrices. Mat. 

Sb. (N.S.), 72(114):4 (1967), 507-536. 
[19] S. Negahban and M. J. Wainwright. Estimation of (near) low-rank matrices with noise and high- 
dimensional scaling. [arXiv:0912.5 100vl (2009) 
[20] A. Rohde, A.B. Tsybakov. Estimation of High-Dimensional Low-Rank Matrices. Ann. Statist. Volume 

39, Number 2 (2011), 887-930. 
[21] M. Rudelson, R. Vershynin. Non-asymptotic theory of random matrices: extreme singular values. 

|arXiv: 1003.2990^ 2 (2010) 

[22] M. Yuan, A. Ekici, Z. Lu and R. Monteiro. Dimension Reduction and Coefficient Estimation in 
Multivariate Linear Regression. Journal of the Royal Statistical Society, Series B, 69 (2007), 329-346. 

Annex A. Monte Carlo evaluation of S qxn (r) 

SMonteCarlo <- function(q,n,Nsim=200){ 
Sk <- array (0 , c(Nsim,min(q,n) ) ) 
for (is in l:Nsim) { 

s <- svd (matrix (rnorm (q*n) ,nrow=q,ncol=n) ,nu=0,nv=0)$d 
Sk [is ,] <-sqrt (cumsum(s**2) ) 

> 

return(apply (Sk, 2 ,mean) ) 



Annex B. Marchenko- Pastur approximation of S qxn (r) 

SMarchenkoPastur <-f unction (q,n, eps= 10** (-9) ){ 
beta <- min (n , q) /max (n , q) 
alpha <- (1 :min(n,q) )/min(n,q) 
s<-rep (0 , min (n , q) ) 
f <- function(x){ 

return (sqrt ( (x- (1-sqrt (beta) ) "2) * ( (1+sqrt (beta) ) ~2-x) ) / (2*pi*beta*x) ) > 
xf <- function(x){ 

return(sqrt ( (x- (1-sqrt (beta) ) "2) * ( (1+sqrt (beta) ) ~2-x) ) / (2*pi*beta) ) } 
for (a in 1 : length (alpha) ) { 

m <- (1-sqrt (beta) )~2 
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M <- (1+sqrt (beta) )~2 
while ((M-m)>eps) { 

if (integrate (f, (m+M)/2, (1+sqrt (beta) ) ~2) $value<alpha [a] ) M<-(m+M)/2 

else m<-(m+M)/2 

} 

s[a] <- integrate (xf , (m+M)/2, (1+sqrt (beta)) "2) $value 

> 

return (sqrt (s*n*q) ) 

} 

CMAP, UMR CNRS 7641, Ecole Polytechnique, Route de Saclay, 91128 Palaiseau Cedex, 
France 

E-mail address: christophe.giraud@polytechnique.edu 
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Figure 3. Experiment 1. Left to right : p = 0.1, 0.5, 0.9. Top : box- 



plots of the ratio Q for KF[K = 2], RSCIfJC = 2], KF[K = CV] and 
RSCI[X = CV}. Bottom : mean estimated rank E[f] for each estimator 
and each value of b. 
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Figure 4. Experiment 2. Left to right : p = 0.1, 0.5, 0.9. Top : box- 



plots of the ratio ([14]) truncated at 10 for KF[iT = 2], ~KF[K = CV] and 
RSC[A = CV]. Bottom : mean estimated rank E[f] for each estimator and 
each value of b. 



